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We extend the nested sampling algorithm to simulate materials under periodic boundary and 
constant pressure conditions, and show how it can be used to determine the complete equilibrium 
phase diagram, for a given potential energy function, efficiently and in a highly automated fashion. 

The only inputs required are the composition and the desired pressure and temperature ranges, 
in particular, solid-solid phase transitions are recovered without any a priori knowledge about the 
structure of solid phases. We benchmark and showcase the algorithm on the periodic Lennard-Jones 
system, aluminium and NiTi. 


I. INTRODUCTION 

Phase diagrams of materials describe the regions of sta¬ 
bility and equilibria of structurally distinct phases and 
are fundamental in both materials science and industry. 
In order to augment experiments, computer simulations 
and theoretical calculations are often used to provide ref¬ 
erence data and describe phase transitions. A plethora of 
methods exist to determine individual phase boundaries, 
including Gibbs ensemble Monte Carlo [1], Gibbs-Duhem 
integration [2], thermodynamic integration and even di¬ 
rect molecular dynamics simulations of coexistence. Each 
of these algorithms requires the user to specify at least 
the identity and approximate location of the phase tran¬ 
sition under investigation. Moreover, in the case of the 
solid phases, where much of the interest lies, advance 
knowledge of the crystal structure of each phase is re¬ 
quired. Calculating an entire phase diagram by combin¬ 
ing the results of such methods therefore demands a high 
degree of prior knowledge of the result. This in turn 
poses a barrier to the discovery of unexpected phases 
and phase transitions. Furthermore, such algorithms re¬ 
quire specific expertise and separate setup for each type 
of phase transition. 

In this paper we introduce a single algorithm, based 
on nested sampling (NS) [3, 4], that enables the efficient 
calculation of complete pressure-temperature phase di¬ 
agrams, including the solid region. This algorithm re¬ 
quires no prior knowledge of the phase diagram, and 
takes only the potential energy function together with 
the desired pressure and temperature ranges as inputs. 
Moreover, the direct output of the simulation is the par¬ 
tition function as an explicit function of its natural vari¬ 
ables, so calculating thermodynamic observables, such as 
the heat capacity, is straightforward. 

Nested sampling systematically explores the entire po¬ 
tential energy landscape, and in this way is related to 
parallel tempering (also known as replica exchange) [5, 6] 
and Wang-Landau sampling [7]. However, those algo¬ 
rithms encounter a particular convergence problem at 
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first order phase transitions because the probability dis¬ 
tributions (parametrised in terms of temperature in case 
of parallel tempering or energy in case of Wang-Landau) 
on the two sides of the phase transition have very little 
overlap [8]. This results in poor equilibration between 
the distributions on either side of the phase transition 
and large errors (both random and systematic) in the 
predicted locations of phase transitions. 

The NS algorithm was designed to solve this problem. 
It constructs a sequence of decreasing potential energy 
levels, {Ei}, each of which bounds from above a volume 
of configuration space yy, with the property that \i is ap¬ 
proximately a constant factor smaller than the volume, 
Xi- 1 , corresponding to the level above. Each volume is 
sampled uniformly, and therefore each distribution will 
have an approximately constant fractional overlap with 
the one immediately before and after, ensuring fast con¬ 
vergence of the sampling and allowing an accurate eval¬ 
uation of phase space integrals. In particular, the energy 
levels near the phase transition, where phase volumes 
change rapidly, will be very narrowly spaced. The se¬ 
quence of energy levels comprise a discretisation of the 
cumulative density of states x{E ), which allows the eval¬ 
uation of the partition function at arbitrary tempera¬ 
tures, 

Z(N,V,I3) = T (^) 3W2 / dEx'We-^ (1) 

^Z m (N,/3)J2(x l -i~Xi)e- f)E ' (2) 

i 

where N is the number of particles of mass in, V is the 
volume, (3 is the inverse temperature, h is Planck’s con¬ 
stant, the density of states x! is the derivative of %, and 
we labelled the factor resulting from the momentum in¬ 
tegral as Z m . The total phase space volume is xo = V N 
corresponding to the ideal gas limit. Note that the se¬ 
quence of energies {Ei} and configuration space volumes 
{Xi} are independent of temperature, so the partition 
function can be evaluated a posteriori at any tempera¬ 
ture by changing /3 in (2). 

The basic NS algorithm is as follows. We initialise by 
generating a pool of K uniformly random configurations 
and iterate the following loop starting at i = 1. 
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1. Record the energy of the sample with the highest 
energy as Ei , and use it as the new energy limit, 
Aliimit 4— Ei. The corresponding phase space vol¬ 
ume is Xi « Xo[K/(K + 1)]\ 

2. Remove the sample with energy Ei from the pool 
and generate a new configuration uniformly ran¬ 
dom in the configuration space, subject to the con¬ 
straint that its energy is less than -Ejimit- One way 
to do this is to clone a randomly chosen existing 
configuration and make it undergo a random walk 
of L steps, subject only to the energy limit con¬ 
straint. 

3. Let i 4— i + 1, and return to step 1. 

At each iteration, the pool of K samples are uni¬ 
formly distributed in configuration space with energy 
E < -Eumit- The finite sample size leads to a statisti¬ 
cal error in logXi, and also in the computed observables, 
that is asymptotically proportional to 1 /\/K, so any de¬ 
sired accuracy can be achieved by increasing K. Note 
that for any given AT, the sequence of energies and phase 
volumes converge exponentially fast (the number of iter¬ 
ations required to obtain results shown below never ex¬ 
ceeded 2000 • AT), and increasing K necessitates a new 
simulation from scratch. 

Since its inception NS has been used successfully for 
Bayesian model selection in astrophysics [9], and also to 
investigate the potential energy landscapes of atomistic 
systems ranging from clusters to proteins [10-18]. 

The structure of this paper is as follows. In section II 
we modify the NS algorithm to enable its application 
at constant isotropic pressure with fully flexible periodic 
boundary conditions [19] where the periodic simulation 
cell is allowed to change shape. In sections III and IV 
we show that this development enables the determina¬ 
tion of pressure-temperature phase diagrams of materi¬ 
als directly from the potential energy function without 
recourse to any other a priori knowledge. In particular, 
in section IV we calculate phase diagrams for aluminium 
and NiTi. Finally in section V we conclude this paper, 
discussing some consequences of the capability to calcu¬ 
late entire phase diagrams with a single method and in 
a highly automated fashion. 

II. NESTED SAMPLING WITH FULLY 
FLEXIBLE PERIODIC BOUNDARY 
CONDITIONS AT CONSTANT PRESSURE 

Nested sampling produces new samples by cloning an 
existing sample and then evolving the clone using a 
Markov chain Monte Carlo (MCMC) random walk [20]. 
Although one could work in the NVT ensemble and use 
equations 1 and 2, that would be very inefficient. MCMC 
simulations performed at fixed pressure require just a 
fraction of the computational expense as equivalent cal¬ 
culations performed at fixed volume. There are two rea¬ 
sons for this. 


First, allowing the system to change volume by dilat¬ 
ing or contracting expedites the cooperative freeing of 
jammed atoms. In contrast, at fixed volume, atoms that 
have become jammed are only freed by the coincidental 
movement of all atoms to separate them. Consequently 
MCMC simulations at fixed pressure explore configura¬ 
tion space far more rapidly than simulations at fixed vol¬ 
ume. 

The second reason arises from the thermodynamic be¬ 
haviour of systems at a first order phase transition. At 
a phase transition under constant volume conditions the 
two phases coexist and an interface forms between them. 
Such interfaces are large on the atomic scale [21] and the 
behaviour of atoms at an interface is not representative 
of the behaviour of atoms in the equilibrium phases. As 
a result the interface introduces a systematic error that 
is only overcome by simulating very large numbers of 
atoms. 

Such interfaces also occur under constant pressure con¬ 
ditions in the infinite system size limit. The contribution 
to the Gibbs Free Energy from an interface is propor¬ 
tional to 7 N 5 , where 7 is the interfacial tension. In 
contrast, the Gibbs Free Energies of each of the pure 
phases are extensive (proportional to N). Therefore the 
Gibbs Free Energy cost of the interface is negligible for 
thermodynamic systems. Conversely, for the relatively 
small system sizes amenable to density of states calcula¬ 
tion methods such as nested sampling, the Gibbs Free En¬ 
ergy cost of the interface is appreciable, provided 7 is not 
close to zero. Consequently, at a constant pressure phase 
transition between phases with identical atomic compo¬ 
sitions, configurations containing an interface have neg¬ 
ligible statistical weight in such simulations, and a dis¬ 
continuous transition is observed from one equilibrium 
phase to the other. This enables the accurate simula¬ 
tion of phase transitions using much smaller numbers of 
atoms. 

Using small numbers of atoms to simulate a phase tran¬ 
sition naturally introduces new finite size errors. In par¬ 
ticular, for a fixed number of atoms, it is not possible to 
represent all crystal structures in a simulation cell of fixed 
shape. This representational bias is removed by the use 
of fully flexible periodic boundary conditions [19], which 
allow the simulation cell to deform smoothly and thus 
take any shape. However, using fully flexible periodic 
boundary conditions allows the formation of very thin 
simulation cells containing unphysical quasi one and two 
dimensional configurations, characterised by interacting 
periodic images. In subsection IIA we describe a rigorous 
solution to this new finite size problem. Later, in sub¬ 
section IIB we describe the calculation of the constant 
pressure partition function and heat capacity, both as 
explicit functions of temperature, using nested sampling. 
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A. Constraint on the simulation cell to exclude 
unphysical quasi one and two dimensional 
configurations 


The partition function at fixed isotropic pressure p 
with fully flexible periodic boundary conditions [19] is 


A (N,p,/3) = Z m /3p J dh 0 S (det h 0 - 1 ) x 


/ dVV N / dse- /3H(s ’ h °’ v> ). 

/o o,i) 3JV 


( 3 ) 


Here H (s, h 0 , V,p) = E (s, h 0 , V) + pV, h is the 3x3 
matrix of lattice vectors relating the Cartesian positions 
of the atoms r to the fractional coordinates s via r = hs, 
V = deth is the volume, and ho = hV ~ 1,/3 is the image 
of the unit cell normalised to unit volume. 

The partition function (3) corresponds to integration 
over all nine elements of the matrix ho, and the (5-function 
restricts the integration to matrices satisfying det h 0 = 1 . 
This partition function is formally correct in the thermo¬ 
dynamic limit [19, 22]. However, finite systems in this 
description can adopt configurations for which the sim¬ 
ulation cell becomes very thin. In this case, periodic 
boundary conditions give rise to a quasi one or two di¬ 
mensional system. The prevalence of such configurations 
leads to a poor approximation of the three dimensional 
atomic system due to excessively large finite size effects. 
We exclude such thin configurations by changing the lim¬ 
its for integration over elements of ho, so that the perpen¬ 
dicular distances between opposite faces of the simulation 
cell h 0 are greater than some “minimum cell depth” value 
do- 

The perpendicular distance between faces of the unit 
cell h made by lattice vectors h^ and h (j -* is given by 


_i_ det h 

h(k) = jhW xhW|' 


( 4 ) 


The cell depth D( ho), which measures how “thin” the cell 
has become, is defined as the minimum value of d ^ (k) , for 
the cell at normalised (unit) volume ho. 


L>(h 0 ) = min 

i= 1,2,3 V K J 


( 5 ) 


Thus we integrate over elements of h 0 such that 

D (h 0 ) > do- ( 6 ) 

The minimum cell depth do is a real number on the in¬ 
terval [ 0 , 1 ] where do = 1 restricts the simulation cell to a 
cube. Smaller values of do are accordingly less restrictive 
on the shape of the simulation cell, and do = 0 corre¬ 
sponds to no restrictions on the simulation cell. 

Incorporating this change of integration limits into the 
partition function (3) yields a new partition function 


A (N,p, j3, do) = Z m /3p / dh 0 S (det h 0 - 1) x 
J D(ho)>do 

rdVV N [ dse-W sMp) - 

Jo J( o,i) 3JV 


( 7 ) 


In the thermodynamic limit (7) is equal to (3) up to 
a factor which depends only on do- The two partition 
functions are equal if and only if do = 0 . 

In tests with 64 atoms we verified that the heat capac¬ 
ity curves were independent of do at values of 0.65, 0.7 
and 0.8, in Lennard-Jonesium and aluminium. The win¬ 
dow of independence from do grows wider as the number 
of particles is increased. For larger numbers of atoms, 
there are more ways to arrange those atoms into a given 
crystal structure, including in simulation cells that are 
closer to a cube. Similarly, unphysical correlations are 
introduced when the absolute number of atoms between 
faces of the cell becomes too small, and therefore larger 
simulations can tolerate “thinner” simulation cells ho. 
The nickel-titanium calculations were performed with 
d 0 = 0.7. 


B. Partition function and thermodynamic variables 

The partition function we seek to calculate is given in 
equation (7). Above some sufficiently large volume Vo, we 
approximate the system as an ideal gas, neglecting inter¬ 
atomic interactions, which corresponds to the condition 
E(s,h 0 ,V) <C pV. In this approximation the volume 
integral in (7) is the sum of two parts 


A(N,p,fi, d 0 ) « Z m /3p 


Ans(^i P, P, Vo, do) 


+ f dh 0 <5 (det h 0 - 1) f dVV N [ dse~ 0pV 
JD(h 0 )>d 0 Jv 0 J( 0,1) 3JV 


( 8 ) 


where 


A NS (N,p,l3, Vo, d 0 ) = / dh 0 S (det h 0 - 1) x 
J D(ho)>do 
rVo r 

/ dVV N / dse -0{E(s,h o ,V)+pV] 

Jo j( 0,1) 3JV 

(9) 

We calculate A N $ using nested sampling. Calculations 
are performed at fixed pressure to generate a sequence 
of enthalpies, Hi, where H = E (s, V, h 0 ) + pV. The NS 
approximation for Ans> is 


A m (N,p,/3,Vo,d 0 ) «E(x»-i — Xi) e l3Hi 

ZL (10) 

« E A Xi e~ pHi 

i=l 

where \i ~ Xo [k+i ) > Xo = and A Xt . ~ X i- i - 

X i- We use single atom Monte Carlo (MC) moves in 
fractional coordinates with the amplitude updated every 
^ iterations to maintain a good acceptance rate. Uni¬ 
form sampling of lattice shape matrices ho subject to 
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equation ( 6 ) was achieved by independent shearing and 
stretching moves which do not change the volume. The 
ratios of atom, volume, shear and stretch moves were 
IV:10:1:1. Further details of the MC moves and paral¬ 
lelisation scheme are given in the Supplemental Material 
(SM) [23]. 

We show in appendix A that volumes greater than 
Vo make a negligible contribution to the partition func¬ 
tion ( 8 ), provided k B T <C pVo. In this case we have 

/ \ 3iV/2 

A (N,p,P,do) w ^ A NS (IV,p,/3, V 0 ,d 0 ) 

(U) 

where we have expanded Z m . One can always assert the 
condition k B T pV o, and in practice it is easy to find 
values of Vo suitable for physically relevant conditions. 
We found Vo = 10 7 N A 3 to be suitable for all condi¬ 
tions considered in this paper. From (11) we obtain the 
expected enthalpy 


(H)=~ 


dlogA(N,p,(3,d 0 ) 

dp 


/31V \ 1 

— I 2 1 ) ^ “t“ (-^configurations/ 


( 12 ) 


and the heat capacity at constant pressure 

= -*■>' r aJ w 

+ k B P ((-^configurations) — (-^configurations) 


(13) 

(14) 


where 


(-^configurations) 


El°r H < e ~ pHt 

Ei=“ Ax* e~P H ' ’ 


E-"i x A Xi Hf e~P Hi 

configurations ) ~ ^i max ^ 


(15) 


(-^configurations ) 


This form (15) naturally does not depend on the contri¬ 
bution made by the low density configurations omitted 
from the NS calculation, or explicitly on the value of do. 
We used equations (15) when calculating the heat capac¬ 
ities presented in this paper. 


III. CALCULATING PHASE DIAGRAMS 

In this section we describe a method for calculating the 
phase diagram of a material from the output of nested 
sampling. We then benchmark the performance of nested 
sampling on the periodic Lennard-Jones system, and find 
nested sampling to be orders of magnitude more efficient 
than Parallel Tempering (PT) for resolving the melting 
and evaporation transitions. 
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FIG. 1. Demonstration of how NS can be used to calculate 
phase diagrams, using the case of the periodic Lennard-Jones 
model. NS calculations are performed at a series of pressures 
and phase transitions are located by peaks of the heat capacity 
curves (blue). The red lines show values from the literature for 
the melting (solid) [24], boiling (dashed) [25] and sublimation 
(dotted) [26] curves. 


Given the partition function (11), phase transitions can 
be easily located by finding the peaks of response func¬ 
tions such as the heat capacity (14). By performing sep¬ 
arate NS simulations at a number of pressures and com¬ 
bining the pressure and temperature values correspond¬ 
ing to the heat capacity peaks one can straightforwardly 
construct the entire phase diagram including all thermo¬ 
dynamically stable phases. This process is illustrated in 
Figure 1. 

In Figure 2 we compare the performance of NS to that 
of PT for calculating the melting and evaporation tran¬ 
sitions. NS provides a reasonable estimate of the melting 
and boiling points using only ~ 10 8 energy evaluations, 
while parallel tempering needs many orders of magni¬ 
tude more computational effort than NS to find the evap¬ 
oration transition and almost two orders of magnitude 
more computational effort to find the melting transition. 
(A similar increase in computational efficiency compared 
with parallel tempering was found for LJ clusters [10] 
and hard spheres [11, 27].) 

Finally, in Figure 3 we show the phase diagram for 
64 particles of Lennard Jonesium as calculated using NS 
with K = 640, L = 1.6 x 10 5 . Comparison with the 
literature phase diagrams for ~ 500 particles confirms 
excellent agreement with the literature values for the 
evaporation transition [25] and also the solid-liquid and 
high pressure solid-vapour transitions [24]. Below the 
triple point, we observe slower convergence with respect 
to L towards literature values of the sublimation transi¬ 
tion [26]. We also find the beginning of the Widom-line: 
the shallow line of heat capacity maxima that extends 
into the supercritical region. The Widom-line and our 
method for estimating the critical point are described in 
the SM [23]. 
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FIG. 2. Performance comparison of NS and PT. Sixty-four 
Lennard-Jones particles were simulated at a pressure of 0.027 
(Lennard-Jones units). Both NS and PT simulations were 
initialised from the vapour phase. PT was performed using 
128 equispaced temperature values in the range [0.4,1.4], The 
left panel shows the estimated transition temperatures as a 
function of computational cost while the right panel shows the 
mean enthalpy as a function of temperature corresponding to 
three selected values of the cost. 



-[24]: 108 particles -[26] 

- [24]: 256 particles —i— NS: phase transitions 
- [25] -+- NS: Widom-line 

FIG. 3. Phase diagram for N = 64 Lennard-Jones particles as 
calculated using NS, with comparison to the literature (N « 
500) phase diagram, as described in the text. 


IV. RESULTS 
A. Aluminium 

In this section we apply the new algorithm to several 
empirical models of aluminium in order to demonstrate 
the capability of nested sampling to find solid-solid phase 
transitions without any prior knowledge of the crystal 
structures or even the existence of multiple stable phases. 
Furthermore, although the particular off-the-shelf mod¬ 
els we use here do not reproduce the experimentally de¬ 
termined phase diagram of the material everywhere, the 
fact that nested sampling allows a direct calculation of 
the entire phase diagram means that in the future one 
could automate the optimisation of potentials to match 


the experimental phase diagram. 

As one of the most commonly used metals, the thermo¬ 
dynamic properties of aluminium have been extensively 
studied. The melting line of aluminium has been mea¬ 
sured up to 125 GPa [28-31], with good agreement be¬ 
tween the different experimental techniques. Theoretical 
calculations have also been performed using embedded- 
atom type potentials [37-44] and ab initio methods [45- 
47], the latter providing melting temperatures up to 350 
GPa [48]. At ambient conditions aluminium crystallises 
in the face-centred-cubic (fee) structure, but a phase 
transition to the hexagonal-close-packed (hep) structure 
at 217 GPa has been revealed by X-ray diffraction exper¬ 
iments [49] and the body-centred-cubic (bcc) phase has 
been also produced in laser-induced microexplosions [50]. 
The critical points of most metals are not amenable to 
conventional experimental study and thus estimation of 
their properties is usually based upon empirical relation¬ 
ships between the critical temperature and other mea¬ 
sured thermodynamic properties. In the case of alu¬ 
minium these result in predictions in a wide temperature 
and pressure range [32-35]. 

We chose four widely used models all based on the 
embedded-atom method (EAM): (1) the model devel¬ 
oped by Liu et al. [43] (LEA-EAM), which is an im¬ 
proved version of the original potential of Ercolessi and 
Adams [42], (2) the model developed by Mishin et al. [44] 
using both experimental and ab initio data (Mishin- 
EAM), (3) the EAM of Mei and Davenport [40] (MD- 
EAM) and (4) the recently modified version of the MD- 
EAM, reparametrised by Jasper et al. to accurately re¬ 
produce the DFT energies for Al clusters and nanoparti¬ 
cles of various sizes (NPB-EAM) [51]. 

The phase diagrams for all four models based on NS 
simulations with 64 particles are shown in Figure 4. The 
resulting critical parameters vary over a wide range for 
the different models. Above the critical point we observe 
the Widom-line, indicated by those points not linked by 
a solid line. Heat capacity maxima corresponding to the 
Widom-line become broader away from the critical point, 
as indicated by the larger error bars. The Widom-line 
and our method for estimating the critical point are de¬ 
scribed in the SM [23]. 

The melting lines are in a good agreement with the 
available experimental data up to the pressure value 
p ss 25 GPa. Above that pressure, the melting curves 
diverge from the experimental results, except for the MD- 
EAM potential, which reproduces the melting curve re¬ 
markably well. 

At higher pressures small peaks appear on the heat 
capacity curves below the melting temperature for all 
models indicating solid-solid phase transitions (see ap¬ 
pendix B). We post-processed the samples from the NS 
simulations. As expected, the fee structure was found to 
be stable at low pressures in all four models. However, 
the models differ markedly in their predictions at high 
pressures. The only commonality between the predicted 
high pressure solid phase diagrams is that the maximum 
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Pressure (GPa) 

FIG. 4. Phase diagrams corresponding to four EAM models of aluminium. Red symbols show the NS results, the error bars 
are calculated as the width at half maximum of the peaks on the heat capacity curves. On the boiling line points are connected 
by a solid line up to the critical point. (The method we used to estimate the critical point is described in the SM [23].) Black 
symbols show experimental melting points measured with Bridgman cells [28], with Diamond anvil cells (DAC (a) [29] and 
(b) [30]) and shock waves (SW) [31]. Different square symbols show estimates of the critical point from experiments, (a) [32], 
(b) [33], (c) [34] and (d) [35]. For NPB-EAM and MD-EAM large black squares show the critical point and smaller black 
squares show the evaporation temperatures, all calculated using Gibbs ensemble Monte Carlo [36]. At pressures below the 
critical point, NS parameters K = 800 and L = 3000 were used (the total number of energy evaluations was 3 x 10 9 for each 
pressure), while runs at pressures where solid-solid transitions are present required K = 3200 and L = 15000 (total number of 
energy evaluations were 4 x 10 10 ). 


predicted stable pressure for the fee structure is far too 
low, both in comparison with experiment and density 
functional theory [49, 52, 53]. 


B. NiTi 

Finally, in order to demonstrate that NS is applicable 
to more complex problems, we show results for a mate¬ 
rial of current scientific interest, the NiTi shape mem¬ 
ory alloy [54, 55]. The shape memory effect relies on 
the structural phase transition from the high temper¬ 


ature austenitic phase to the low temperature marten¬ 
sitic phase [56]. Studying this transition is particularly 
challenging with traditional free energy methods because 
the austenitic phase does not correspond to a local mini¬ 
mum of the potential energy surface. Figure 5 shows the 
pressure-temperature-composition phase diagram corre¬ 
sponding to a recent EAM model [57, 58] as computed 
with NS. The phase transition temperatures are within 
50 K of the experimental values and reproduce the trend 
with compositional change. We predict a decreasing tran¬ 
sition temperature with increasing pressure. It is notable 
that this potential successfully reproduces the marten- 
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FIG. 5. NiTi martensitic phase transition as a function of Ni 
content (at 0.66 GPa) and pressure (at 1:1 composition). The 
simulation cell contained 64 atoms in the cases of the 50% and 
51.6% Ni compositions and 108 atoms in the case of 50.8% Ni 
content. NS parameters were K = 1920, L = 10’, each data 
point used 10 10 energy evaluations, and Ni-Ti swap moves 
were also included in the MC. Experimental results are taken 
from [59]. 


we suggest NS is eminently suitable for validating mate¬ 
rials models, and in the future could even play a role in 
the automatic optimisation of empirical models. 
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Appendix A IDEAL GAS CONTRIBUTION TO 
THE PARTITION FUNCTION 


sitic transition temperature, despite the fact that the 
minimum enthalpy structure for the potential is differ¬ 
ent to the structure observed both experimentally and in 
DFT: here the lowest enthalpy structure (which we label 
B19X) is orthorhombic (see the SM [23] for a descrip¬ 
tion of the low enthalpy structures we identified). Thus 
it appears that the austenite-martensite transition tem¬ 
perature is not sensitive to the detailed geometry and 
ordering of the lowest enthalpy structures. Such empir¬ 
ical potentials can therefore be useful tools for studying 
this transition in the future. 


V. CONCLUSION AND OUTLOOK 

In summary, we have extended the nested sampling al¬ 
gorithm to allow simulations using fully flexible periodic 
boundary conditions at fixed pressure and demonstrated 
how it can be used to determine pressure-temperature- 
composition phase diagrams. In contrast to existing 
methods for comparing specific phases, NS explores the 
entire configuration space without requiring any prior 
knowledge about the structures of different solid phases 
with the only necessary input being the composition 
and the desired pressure and temperature ranges. This 
makes it the method of choice for exploring the pressure- 
temperature-composition space, which is the next un¬ 
explored realm naturally following much recent work in 
crystal structure exploration at zero temperature. Since 
the algorithm is run independently for different pressures 
and compositions, and also has excellent parallel scaling 
up to a number of processors equal to the number of si¬ 
multaneous samples, it might even be possible to run it 
on ab initio models on exascale computers. Furthermore, 


In this appendix we show that the ideal gas contri¬ 
bution to the partition function (8) asymptotically ap¬ 
proaches zero for any positive minimum cell depth do, in 
the limit kBT/pVo —> 0. 

The ideal gas contribution to the partition function (8) 
is 

[ dh 0 <5 (det h 0 — 1) f dVV N [ dse~ 0pV 

JD(h 0 )>d 0 J Vo 0,1) 3JV 

(16) 

We begin by noting that the exponential term does not 
depend on E( s, h 0 ,U), and therefore /( 0 1 \ 3 jv ds = 1. 
Thus we have 

[ dM(detho-l) f dVV N [ dse -ft,v 

JD(h. 0 )>d 0 JVo J (0,1) 3JV 

= [ dh 0 <5 (det h 0 — 1 ) [ dVV N e~ l3pV . 

JD(h 0 )>d 0 Jvo 

(17) 

The integral over volume V evaluates to 

r°° i 

J dW N e~P pV = + 1, PpV 0 ) (18) 

where T(N + 1 ,/3pVo) is the upper incomplete gamma 
function. 

Finally we define the function A (do) to be equal to the 
integral over ho 

A(d 0 )= [ dh 0 <5 (det h 0 — 1). (19) 

J D(h 0 )>d 0 

The function A (do) is finite for any positive value of d 0 , 
A( 1) = 0 and A (do) diverges in the limit do —> 0. In 











the orthorhombic case, where all angles of the simulation 
cell are equal to f, A (do) = | (logdo) 2 , with A = 1 
at do « 0.62. However at any positive value of d 0 the 
contribution to the partition function (8) due to volumes 
greater than Vo goes to zero in the limit k^T/pYo —> 0 
because T(N + 1, /3pV o) — > 0 in the same limit. 


Appendix B IDENTIFYING SOLID-SOLID 
PHASE TRANSITIONS 


The locations of phase transitions are determined 
solely by looking at the peaks in the heat capacity. Next, 
we inspect the system at temperatures either side of the 
phase transition. Specific phases can be identified in the 
following way. If no appropriate order parameter is to 
hand, then one picks a number of random configurations 
from the output of nested sampling, chosen according to 
their thermal weights A Xi e ~^ Hi i and inspects them by 
eye. If an appropriate order parameter is known, one can 
compute the free energy landscape for that order parame¬ 
ter. Here one proceeds by binning the weights A Xi e ~^ Hi 
of all configurations, according to the order parameter, 


to create a partial sum Aj = )U Ay.;e ,3Hi for each bin 
j. The free energy for each bin can then be computed as 
F j = l°g( A j) + log log ■ In fact, 

simply calculating the expected enthalpy at the phase 
transition, and then examining the order parameter val¬ 
ues for output configurations around that enthalpy is of¬ 
ten sufficient to identify the crystal structures. 

An example of the latter approach is shown in Figure 6 
for the Mishin-EAM potential, which compares the en¬ 
thalpies and Qe bond order parameter values for nested 
sampling output configurations at three different pres¬ 
sures. At p = 25.0 GPa no phase transition occurs, and 
only fee configurations are present. At p = 34.9 GPa 
a first order phase transition occurs at the average en¬ 
thalpy marked by the vertical dashed line. At that en¬ 
thalpy there is a clear transition between two basins, from 
a first basin that corresponds to the bcc structure, to a 
second that corresponds to the hep structure. Finally, at 
p = 37.5 GPa no phase transition occurs and so there is 
no peak in the heat capacity. At this pressure the bcc 
structure is stable at all temperatures below the melting 
point. Nevertheless, the hep structure is clearly visible 
as a metastable structure. 
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thalpies, thus from right to left in each plot. Horizontal dotted 
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pected enthalpy at the solid-solid phase transition, which was 
located by inspecting the heat capacity and observing a peak. 
These results are discussed in the text. 
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I. DETAILED MONTE CARLO SCHEME 


To understand the Monte Carlo scheme for NS in the fixed pressure ensemble with a fully flexible cell, it is helpful 
to first consider the Monte Carlo scheme for the fixed volume ensemble, with a fixed cell. In that case the partition 
function is given by 


Z(N,Va3) = ^ (^fV/ dse-WV (SI) 

As described in previous papers [S1-S3], each iteration begins with K samples distributed uniformly in s, with 
E < -Eiimit- That is to say, those K samples are drawn from the distribution 


Prob (s|£i imit ) 


x(-Eiimit)’ ^ ( S ’ 1^) ^ -®limit 

0, Elsewhere. 


(S2) 


For the fixed pressure ensemble, we must explore volume and cell matrices ho, as well as the fractional coordinates 
s. In particular, we calculate the integral (9) with NS. 

[' /' C r 

Ans(N,p,/3, Vo,do) = / dh 0 <5 (det h 0 — 1) x / dVV N / d se - 0 [ E ( s ’ h o,v)+pV] (9, revisited) 

JD(h 0 )>do Jo J(0,1) 3N 

In (9) the volume has a weight V N with V < Vo, and h 0 is constrained to a surface of determinant 1, with D (h 0 ) > do- 
Therefore we seek to generate new samples according to the distribution 


Prob (s,h 0 ,V|Hi im it) 

H (s,h 0 ,V) < V<V 0 , 

= \ D (ho) > do 

[d, Elsewhere. 


(S3) 


As described in the paper, this can be achieved by Markov chain Monte Carlo (MCMC) sampling. MCMC walks 
were random sequences of the following Monte Carlo (MC) moves, which were chosen according to the probabilities 
(ratios) given in the paper, unless specifically stated in this Supplemental Material. 

Atom moves: Single atom MC moves are performed as follows. 

1. Select a single atom at random. This atom has fractional coordinates 

2. Displace that atom’s fractional coordinates by a random vector 

Si —>■ Si -(- L S A s (S4) 

Here A s is a 3-tuple with elements drawn independently from a uniform distribution on the interval [—1,1], 
and L s is a constant that controls the step size. The size of L s was updated periodically during the NS 
calculation, as described in the next subsection, “Updating MC step lengths”. 

3. Calculate the enthalpy L/triai = E (s, h 0 , V) + pV for the system after this atom displacement. 

4. Accept the displacement if fftriai < -Hiimit- Otherwise, the old configuration is kept. 

Volume moves: Volume MC steps are as follows. 

1. Propose a volume move from V to V' where 

V' = V + L V A V (S5) 

Here A y is a random number drawn from the uniform distribution on the interval [—1,1], and Ly is a 
constant that controls the step size. The size of Ly was updated periodically during the NS calculation, 
as described in the next subsection, “Updating MC step lengths”. 
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2 . 

3. 

4. 


The displacement is accepted with probability min 1, 


1 N 


. If the displacement is accepted, continue 


to step 3. Otherwise, the old configuration is kept and the move ends. 

If the displacement was accepted in step 2, calculate the enthalpy H tlia \ = E (s, h 0 , V') +pV' for the system 
after this volume displacement. 

If //trial < //limit the displacement is accepted. Otherwise, the old configuration is kept. 


Lattice shear moves: The unit simulation cell ho is sheared, subject to constraints on both the cell depth and 
the enthalpy. The algorithm is given in Algorithm 1. The size of shear steps is controlled by the parameter 
step_l_sh. This parameter was updated periodically during the NS calculation, as described in the next 
subsection, “Updating MC step lengths”. 


Lattice stretch moves: The unit simulation cell ho is stretched at fixed unit volume, subject to constraints on both 
the cell depth and the enthalpy. The algorithm is given in Algorithm 2. The size of stretch steps is controlled 
by the parameter step_l_st. This parameter was updated periodically during the NS calculation, as described 
in the next subsection, “Updating MC step lengths”. 

(Biliary) atom swaps: For the NiTi alloy, we included binary (Ni- Ti) atom swaps as an additional MC proposal. 
The algorithm for this MC proposal is as follows: 


1. Choose a random Ni atom, and a random Ti atom. 

2. Swap the fractional coordinates of those two atoms. 

3. Calculate the enthalpy H tria \ = E (s, h 0 , V) + pV for the system after the swap. 

4. Accept the swap if //trial < //limit- Otherwise, the old configuration is kept. 


A. Updating MC step lengths 

The values of the step size parameters, L s , Ly, step_l_sh and step_l_st, were updated every y iterations (recall 
that K is the number of samples employed in the NS calculation). This corresponds to an expected reduction in the 
enclosed configuration space volume y to a factor of e ~U Each step size parameter was updated by performing a 
short MCMC exploration using only the corresponding type of MC proposal. The step size parameters were updated 
to maintain an acceptance rate in the interval (0.2, 0.3), and as close to 0.25 as possible, for each MC proposal type. 

Since the calculation was performed in parallel, when collecting data about the acceptance rate for each MC proposal 
type, each processor ran MCMC trajectories with separate, random configurations from the current sample set. In 
general, this avoids setting step lengths appropriate to unrepresentative regions of the bounded uniform distribution. 


II. PARALLELISATION SCHEME 

In order to be able to benefit from supercomputing facilities efficiently, the NS algorithm has to be parallelised. 
Fortunately, possible parallelisation schemes come quite naturally from the method, and here we describe a scheme for 
parallelising within each iteration [S4]. While in the serial algorithm, the new clone is decorrelated using a trajectory 
of walk length L in each iteration, in our scheme parallelised over N proc processors, the new clone is decorrelated 
using a trajectory of walk length L/N proc , along with ( N pmc — 1) randomly chosen configurations which are also 
(independently) propagated through L/N pmc steps. Thus, rather than propagating a single configuration for L steps, 
we propagate N pmc configurations for L/N proc steps. 


III. FINITE SIZE EFFECT 

In order to estimate the effect of using only 64 particles in our simulations we performed tests with 32 and 128 
particles as well for both the periodic Lennard-Jones and the LEA-EAM aluminium potential models. The heat 
capacities of these are compared in Figure SI. The Lennard-Jones test was performed at p = 0.027e/<r 3 , below the 
critical point, hence both the evaporation and melting peaks are present. The aluminium test was performed at 
p = 11.6 GPa, above the critical point, thus only the melting peak is present. The figure shows that increasing the 
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Algorithm 1. Lattice shear MC step. Samples the surface det(ho) = 1 subject to the enthalpy and minimum cell depth 
constraints. ranf() is a random number uniform in [0,1]. 


number of atoms, the heat capacity peaks become narrower as expected, and the melting and evaporation transitions 
shift to lower and higher temperatures respectively. Of course, it would be possible to perform a finite size scaling 
analysis from the results of multiple Nested Sampling calculations, performed at different system sizes. However, 
we felt that the phase diagrams of the various aluminium models differed from each other to the extent that it was 
sufficient to give a rough indication of error based on the full-width-at-half-maximum of the peaks. 
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Algorithm 2. Lattice stretch MC step. Samples the surface det(ho) = 1 subject to the enthalpy and minimum cell depth 
constraints. ranf() is a random number uniform in [0,1]. 


IV. LENNARD-JONES POTENTIALS 


The benchmark comparison to parallel tempering was performed using the truncated and shifted Lennard-Jones 
potential [S5]. Other Lennard-Jones calculations were performed using the truncated but not shifted potential with 
the usual mean-field correction [S5] beyond the cut-off radius. All Lennard-Jones calculations employed a cut-off 
radius r c = 3 ct. 
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FIG. SI. Magnitude of the finite size effect in case of the Lennard-Jones potential (left panel) run with parameters K = 672 and 
L = 15120, and LEA-EAM aluminium potential (right panel) with run parameters K = 1608 and L = 10080, as demonstrated 
by the heat capacity curves of systems of 32, 64 and 128 particles. 


V. ALUMINIUM 

A. Simulation details 

The simulation cell contained 64 aluminium atoms. The same ratio of the different types of MC moves was used 
at every pressure. Sampling is robust with respect to the ratio of MC moves of each type, and for these aluminium 
calculations the ratio of MC moves (atom:volume:shear:stretch) was (16:1:1:1), rather than (N:10:l:l). The minimum 
allowed cell depth, do, was set to 0.65. Spline parameters for the LEA-EAM and Mishin-EAM potentials are only 
provided for configurations where no two atoms are closer than a certain, small radius. For separations smaller than 
this radius we made linear extrapolations to the splines. The ultra-close atomic proximity region is only relevant to the 
ultra-high temperature gas phase, and does not contribute to the thermodynamics of aluminium in the temperature 
and pressure ranges considered in this paper. Therefore the precise form of this extrapolation does not affect any of 
the results we show. 


B. Heat capacity curves 

Heat capacity curves calculated with Nested Sampling are shown in Figures S2, S3, S4 and S5. 


C. Critical point and the Widom-line 

To locate the critical point, we draw on the results of Bruce and Wilding [S6]. Specifically we begin from their 
analysis of the probability distribution for density for finite systems, in the pressure-temperature region that would 
contain the critical point in an infinite system. For a finite system at and below the critical point the density 
distribution for temperatures corresponding to maxima of the heat capacity appears as a bimodal distribution. Above 
the critical point (along the Widom-line), the distribution transitions quickly to a unimodal distribution. In the 
context of a phase diagram, our location of the critical point appears quite precise, but in fact our sampling of the 
pressure axis is relatively sparse compared to the width of the scaling region around the critical point, since we are 
inspecting the system across pressures spanning many orders of magnitude. Therefore, between adjacent pressure 
values the density distribution transitions rapidly between being bimodal and unimodal. We chose the lowest pressure 
for which the distribution is clearly unimodal as an upper bound for the critical pressure. The corresponding upper 
bound for the critical temperature was taken to be the temperature at which the heat capacity, Cp, is its local 
maximum at that pressure. 

According to textbook definitions, beyond the critical point a single fluid phase is defined. However, response 
functions such as the heat capacity, thermal expansion coefficient and compressibility continue to exhibit maxima into 
the supercritical region. Lines of these maxima, which spread out from the critical point are called the Widom-lines. 
If one moves away from the critical point, the different Widom-lines rapidly diverge from one another, and the maxima 
themselves become quickly smeared and disappear, in case of the heat capacity at T ss 2.5T C and p ss 10p c [S7]. The 
same tendency can be observed in our Nested Sampling results, see Fig. S2. 
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FIG. S2. Aluminium heat capacity curves at lower pressures up to the end of the Widom-line for the different studied potential 
models. The arrows indicate the peak of the c p curves closest to the critical point, and in order to show the broadened peaks 
above that, the curves in the supercritical region are shown on a different scale (shown on the right hand side). The size of the 
live set and the total number of MC steps during one iteration was K = 804 and L = 3120, respectively. 



FIG. S3. Aluminium heat capacity curves above the critical point for the LEA-EAM potential. The insets show the solid-solid 
phase transition peaks in a larger scale. The size of the live set was K = 3216 and the total number of MC steps during one 
iteration was L = 15120, or if otherwise, it is indicated in the legend. 


By calculating the expected density as a function of temperature, the equation-of-state can be calculated using 
Nested Sampling. These are plotted for the low pressure region together with Gibbs ensemble Monte Carlo simulation 
results from Ref [S8] in Figure S6. The critical parameters for the different aluminium potentials determined by 
Nested Sampling are summarised in Table SI. 


VI. NITI 

Here we discuss the lowest enthalpy structures we identified for the EAM potential used to describe NiTi [S9]. The 
enthalpic and energetic ordering of these structures is the same for all pressures considered in this paper. Therefore 
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FIG. S4. Aluminium heat capacity curves above the critical point for the NPB-EAM potential. The inset shows the solid-solid 
phase transition peaks in a larger scale. The size of the live set was K = 3216 and the total number of MC steps during one 
iteration was L = 15120. 



FIG. S5. Aluminium heat capacity curves above the critical point for the MD-EAM potential. The size of the live set was 
K = 3216 and the total number of MC steps during one iteration was L = 15120. 


we shall restrict our discussion to their energy values. 

The EAM potential gives rise to a number of low energy structures that show some similarity to the experimentally 
observed low temperature structure (B19’) and to other structures found using DFT (namely B19 and BCO). The 
lowest energy structure of the EAM potential is orthorhombic, but with significant displacement of the atoms from 
their high temperature average B2 positions. Table S2 below gives the structural parameters of the three structures 
with the lowest energies. In addition, we found that a multitude of local minima exist between just 10 and 15 
meV/atom above the lowest energy structure. The symmetries of these configurations were identified using the 
findsym package [S10]. 


[SI] L. B. Partay, A. P. Bartok, and G. Csanyi, J. Phys. Chem. B 114, 10502 (2010). 


TABLE SI. Critical temperature and pressure as estimated from the Nested Sampling calculations for the different potential 
models. Numbers in parenthesis show GEMC results from Ref. [S8] for comparison. 



T c (K) 

Pc (GPa) 

LEA-EAM 

6810 ± 170 

0.262 

Mishin-EAM 

5800 ± 140 

0.054 

NPB-EAM 

6100 ± 150 (6299 ±48) 0.083 (0.0896 ± 0.0019) 

MD-EAM 

3340 ± 140 (3381 ± 13) 0.042 (0.044 ± 0.002) 




















FIG. S6. Density-temperature phase diagram of aluminium with the NPB-EAM potential. Red lines show the Nested Sampling 
results at different pressures: 0.000402 GPa, 0.0012 GPa, 0.0208 GPa, 0.0337 GPa, 0.0503 GPa, 0.0679 GPa, 0.0830 GPa, 0.100 
GPa, 0.115 GPa, 0.150 GPa, 0.230 GPa, 0.399 GPa, 0.807 GPa, 1.0 GPa, 2.5 GPa. The curve corresponding to the critical 
pressure is shown by a dashed line. Black and open circles show the liquid and vapour densities respectively, from [S8], with 
a black star showing their estimate of the critical point. The inset shows the minimum value of the derivative of the curves 
obtained by Nested Sampling, with the arrow indicating the pressure value where the plateau between vapour and liquid phases 
diminish at the critical point. 


TABLE S2. Low energy structures of NiTi according to the EAM potential. Lattice vectors are in A, energies are per atom 
and are relative to the B2 phase. 


Symmetry 

a b c 

P 

Wyckoff positions 

Energy (meV) 

Pm3m (B2) 

3.0098 3.0098 3.00980 

90 

Ni: a, Ti: b 

0 

Pmm2 (B19X) 4.02186 4.40909 3.01484 

90 

Ni: b (z=-0.04073), c (z=0.04305) Ti: a (z=-0.4434), d (z=0.4458) 

-53.8 

C2/m (B19) 

4.46697 4.02169 3.00532 

98.1 

Ni: a, Ti: d 

-51.1 

P2i/m (BCO) 2.99360 4.00083 4.88392 

107.8 

Ni: e (x=0.08755, z=-0.3249) Ti: e (x=0.3583, z=0.2167) 

-46.8 
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